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I. INTRODUCTION 

Orderliness, reflected in either scaling properties^ or 
power laws^^l, is encountered in different frameworks in- 
volving social groups. Salient example is Zipf's law^, a 
power law with exponent —2 for the density distribution 
function that is observed in describing urban agglomer- 
ation^ and firm sizes all over the worlcP. This kind 
of "order" has got immense attention in the literature. 
The pertinent regularities have been found in other sce- 
narios as well, ranging from percolation theory and nu- 
clear multi-fragmentatiorP to the abundances of genes in 
various organisms and t issues^, the frequency of words 
in natural language s] 5 * 10 1 , the scientific collaboration net- 
works^, the total number of cites of physics journals^, 
the Internet traffic^ or the Linux packages links^. An- 
other special regularity in the density distribution func- 
tion of the number of votes in the Brazilian elections,^ a 
power law with exponent —1. This kind of distribution is 
the origin of the so-called First-Digit Law, or Benford's 
Law^. Such exponent has also been found in Ref. [18] by 
recourse to an information-theoretic methodology^, for 
both the city-size distribution of the province of Huelva 
(Spain) and the results of the 2008 Spanish General Elec- 
tions. Also log-normal distributions are applied to cities' 
sizes. While a city-size distribution is often associated 
with Zipf's law, this holds only in the upper tail. How- 
ever the entire distribution of cities, not just the largest 
ones, is close to log-normal (see Ref. [17]). Undoubtedly, 
some aspects of human behavior reflect a kind of "uni- 
versality" . 

Common to all these disparate systems is the lack of a 
characteristic size, length or frequency for the observable 
of interest, which makes it scale-invariant. In Ref. [TB] 
we have introduced an information-theoretic technique 
based upon the minimization of Fisher's information 
measure^ (abbreviated as MFI) that allows for the for- 
mulation of a "thermodynamics" for scale-invariant sys- 
tems. The methodology establishes an analogy between 
such systems and physical gases which, in turn, shows 



that the two special power laws mentioned in the pre- 
ceding paragraph lead to a set of relationships formally 
identical to those pertaining to the equilibrium states of a 
scale invariant non-interacting system, the scale-free ideal 
gas (SFIG). The difference between the two distributions 
is thereby attributed to different boundary conditions on 
the SFIG. 

Many social systems, however, are not easily included 
into the two types of description of the prec eding pa ra- 
graph. Inspired in opinion dynamics models^ 5 * 20 * 21 ! we 
described in Ref. [22] a numerical process that reproduces 
the shapes of the empirical city-population distributions. 
The model is based in a competitive cluster growth pro- 
cess inside an scale-free ideal network (a scale-free net- 
work which degree distribution is described as an SFIG). 
The finite size of the network introduces competition, and 
thusly correlations in the cluster sizes. The larger the 
competitiveness, the larger the deviation from the SFIG 
distribution. The equivalence with the workings of phys- 
ical gases is compelling indeed. However, some aspects of 
the concomitant problems defy full understanding, since 
an analytical prescription for the classification of the size- 
distribution of social groups was missing. In order to do 
look for such analytical understanding we will appeal here 
to a variational approach centered on information theory. 
We have recently shown (see Ref. |2"3|) how to deal with 
the well-established MaxEnt techniqu e 25 * 26 ! by including, 
within the associated Lagrangian, information regarding 
equations of motion, i.e., "dynamical" information of a 
kind that goes beyond the customary one, based upon 
expectation values. We demonstrated that proportional 
growth and hyper-exponential growth can both be de- 
scribed in dynamical terms, accurately predicting the 
equilibrium distribution for such systems. This is tan- 
tamount to giving a dynamical interpretation to the 'in- 
formation cost' introduced in Ref. 1271 Our main present 
purpose is to apply what we learned in [53] to empirical 
systems, this demonstrating i) the applicability of Max- 
Ent to collective human behavior, and ii) the potential 
use of an explicit social thermodynamics. 
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This work is is organized as follows: in Sec. II we revisit 
the theoretical MaxEnt approach to be followed here. 
In Sec. Ill we present empirical data for the Spanish 
provinces. The ensuing results clearly show the useful- 
ness of the theoretical method. In Sec. IV we discuss 
some noteworthy features of the results an their possible 
uses in other contexts. A summary is given in Sec. V 
and, finally, some technical aspects are the subject of the 
Appendix. 



II. THEORETICAL APPROACH 
A. Basics 

Following the methodology described in[23l let us define 

i) N as the total population, 

ii) n c as the total number of cities into which people are 

apportioned, 

iii) Xi as the population of the z-th city, 

iv) Xq as the minimum possible population amount 

(which is at least Xq = 1), 

v) p(x) as the relative number of cities with exactly a 

population amount given by x. 

Considering the continuous limit of the distribution p(x), 
the conservation of both N and n c guarantees 



i) / dxp(x) = 1, and 

Jxq 



ii) / dxp{x)x = N/n c . 

We now turn on MaxEnt: in dynamical equilibrium, the 
distribution p(x)dx is the one that maximizes Shannon's 
entropy S subject to the above constraints. One thereby 
concocts the variational Lagrangian 



H = S- /m c - XN, 



(1) 



with fj, and A the corresponding Lagrange multipliers. 
MaxEnt entails SH = 0. 



B. Dynamics 

We showed in Ref. [22] that cities' population expansion 
can be modelled by means of clusters' growth in complex 
networks. Such process (diffusion) in networks generally 
starts i) using a node as a seed, ii) its first neighbors being 
added to the cluster in the first iteration, iii) the neigh- 
bors of those neighbors afterwards, and so on. Mathe- 
matically, if At is the interval of time for each iteration 



and x(t) is the population of the cluster at the time t, we 
can write 



Xi(t) 

Xi(t + At) = Xi(t) + At ^2 k i 



(2) 



where kj is the number of neighbors of the j'-th node 
accepted to the cluster per unit time. According to the 
central limit theorem, we generalize this equation as 

H (t + At) ~ Xi (t) + At (kx, (t) ± a k t/x$)) > ( 3 ) 

where k is the global mean value of accepted neighbors 
per unit time and its standard deviation. For large 
enough values of the population this last term can be 
neglected, and write for the continuous limit 



x l (t) = kxi(t), 



(4) 



where the dot represents the time derivative. This equa- 
tion represent the so-called proportional growth, whose 
generalization to a random process leads to the so-called 
geometrical Brownian motion. We show in Ref. |2"51 how 
to proceed to construct the entropy S for such a scale-free 
dynamics. 



C. MaxEnt for scale invariant systems 

We now consider a system governed by the proportional 
growth of Eq. Q, with k being regarded as reflecting 
a Wiener process. Our equation is linearized by intro- 
duction of the new variable u = \og(x/xo), which leads 
to u = k, and then to a Shannon measure expressed in 
terms of u 



S = 



p(u) \og{p{u))du. 



(5) 



With our two conservation rules written now in u— terms, 
the MaxEnt variational problem becomes 



SH 



S 



/i / dup(u) 
lo 



A / dup(u)e v 
lo 



0, 



(6) 

where we have employed the definition A = xo\. 
The solution is the equilibrium density p{u)du = 
Z~ x exp(— Ae u )du, that one now recasts in terms of the 
observable x getting 



p(x)dx = Z 



e -Ax/x 

L dx. 



(7) 



From the conservation rules we easily obtain for the per- 
tinent constraints the values 



z = r(o,A) 

e~ A N 



(8) 



Ar(0,A) u c xq 
with r(a, z) the so-called incomplete Gamma function. 
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D. The Gamma Scaling Law 

The concomitant cumulative function P(x) reads 

r{0,Ax/x ) 



P{x) = 1 - 



r(o,A) 



(9) 



and the associated rank-distribution (RD) — obtained 
from the inversion of the cumulative one — becomes 



a; = ^r- 1 [o,r(o 1 A)r/n c ] 



(10) 



where r is the (continuous) rank from to n c , and T _1 (z) 
denotes the inverse function of V : T(T^ 1 (z)) = z. We 
derive from here what we call the Gamma Scaling Law 



x' = xA/xq 
r' =rr(0,A)/n c , 

obtaining a "scaled" RD 



(11) 



x> = r-'(o, r i), (12) 

which no longer depends on N, n c , or Xq. 

E. Beyond proportional growth 

We now consider the more general expression 

x(t) = kx q (t), (13) 

where q is a dimensionless parameter. Linearization 
is achieved by means of the q-logarithm of Tsallis' 
statistics^! 



log g (x) = - 1 ; [x > 0; log! (as) =lnat], 
e q (x) = [l + (l- q )x] 1 / {1 - q \ 



(14) 



defining the new variable w = \og q (x / xq) , and so obtain- 
ing w = k. The 'V-Shannon entropy becomes 

/>oo 

S = I p(w)\og(p(w))dw. (15) 



The concomitant MaxEnt problem is of the form 



<)[] =<) S - ii I daiA "' ) — \ I dup(w)e„(w) = 0, 

(16) 

whose solution is p(w)dw = Z 1 exp(— Ae q (w))dw. 
Changing back to the observable x we find 



p(x)dx = Z 



-Ax/xo 

- 1 dx. 

xi 



(17) 



The constraints derive from the conservation rules in the 
usual manner as 



Z=(A/x ) q - 1 T(l-q,A) 
r(2-qr,A) _ N 
Ar(l-<j,A) n c x ' 



(18) 



The current cumulative function turns out to be 

_ r(l-g, As/sp) 
() r(l-«,A) ' [ ' 

and for the associated RD has 

x=^T-i[l-q,T(l-q,A)r/n c ]. (20) 

F. Proportional drift in q-exponential growth 

We now consider a q-equilibrium system subject to a slow 
proportional drift, that may account for the "natural" 
population-growth or fluctuations of proportional nature. 
One is here including a kind of noise that affects the 
underlying dynamical equation via 



X — k\X ~\~ kqX^ j 

or, in terms of u = log(ir), 

u = h + k q e {q ~ 1)u , 



(21) 



(22) 



where k\ characterizes proportional drift and k q hyper- 
exponential growth. Considering k\ also as a Wienner 
process, the system departs equilibrium via proportional 
diffusion. As seen in Ref. UHl the kernel function for that 
kind of diffusion is a log-normal distribution. We face 
the convolution 



p(x; A, xq, q, cr)dx = 
dx 

(A/x )«-ir(l-< z ,A) 



dx' 



e -Ax'/x p-log (x'/x)/2a 



x'i 



X\Jll 



(23) 

where a is the drift-scale parameter. The mean value 
reads 



r(2-g,A) ca 2 /2 = N ^ 
AT(l-q,A) n c x ' 



(24) 



III. EMPIRICAL OBSERVATIONS 

A. Description of the data used in this work 

The raw data are obtained from the Spanish state insti- 
tute INE and cover annually the period 1996-2010 (with 
the exception of 1997). It encompasses up to 8000 munic- 
ipalities (the smallest Spanish administrative unit) dis- 
tributed within 50 provinces (the building blocks of the 
autonomous communities, equivalent of the USA-states). 
The autonomous cities of Ceuta and Melilla are not in- 
cluded. We use provinces and municipalities as the clos- 
est representatives of the ideal of a closed system's funda- 
mental elements. But some arbitrariness remains in the 
data, since i) there are many municipalities that actually 
contain more than one population-center, ii) provinces 
that contains more than a single socioeconomic cluster, 
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and iii) economic regions that extend beyond the fron- 
tiers of a province. Those facts introduce systematic er- 
rors into the data, but we retain enough accuracy for the 
purposes of this work. 

The appropriateness of provinces as the proper scale for 
the analysis of this work is due to a trade-off between 
large enough data sets and locality. Autonomous commu- 
nities as a whole provide a large enough sample, but the 
statistical properties we are interested in change within 
them. Smaller administrative units called comarcas have 
a more local nature, but only contains a few municipali- 
ties and thus introduce large numerical errors. We con- 
sider that choosing the adequate scale is the most im- 
portant ingredient for an study of the kind we attempt 
here. 

The data cover a total Spanish population of 39106917 
inhabitants in 1996, reaching 47254510 in 2010. For 
this last year, the largest municipality (Madrid) covers 
3273049 people and the smallest one just 5 persons (Illan 
de Vacas, Toledo). The total population of each province, 
N, ranges between 6458684 inhabitants (Madrid) and 
95258 (Soria), and the number of municipalities, n c , 
between 371 (Burgos) and 34 (Las Palmas). These 
n c — figures are not large enough to build up an accurate 
density distribution. Accordingly, we have systematically 
worked with RD instead. RD are usually built from the 
vector x = {^j}"=ii where Xi is the population of the 
i-th municipality. We assign rank numbers ranging from 
1 to n c , from the largest (max(x) — > r = 1) to the lowest 
(in population) (min(x) — > r = n c ). In order to compare 
with theoretical, continuous rank-distributions, we have 
found it to be more accurate to assign "middle-point" 
values from r = 0.5 to r = nc — 0.5, instead. 
For our analysis we generate, for each province and for 
each time-period, the pair (ui,iii), with the logarithm of 
the population correlated with the relative population- 
increment iii. They are computed as 



ui = \og[x l {t 2 )x i {ti)\/2 
. _ \0g[Xi{t 2 ) / x l (t 1 )] 



(25) 



to-tl 



where t\ and ti are consecutive years for which data are 
available. 



B. Zeroth order approximation: proportional 
growth 

1. A paradigmatic example: Alicante 

We use as paradigmatic example the province of Alicante 
to illustrate our empirical analysis. The inset in the top 
panel of Fig.[T]depicts u versus ii for all our municipalities 
in the last 15 years. The correlation coefficient between 
the two variables is found to be 0.16, small enough to 
consider, at zeroth order approximation, that the growth 
is of proportional nature, independently of the size, thus 
obeying Eq. Q. 



As seen in the previous section, the equilibrium rank- 
density predicted for proportional growth by MaxEnt fol- 
lows Eq. ([T]) , where A is univocally determined by Eq. Q 
with i) an estimated minimum size xq, ii) the given total 
population N = 1926285 (in 2010) and, iii) the number 
of municipalities n c = 141. N and n c are well deter- 
mined in all cases. In assigning xq, fluctuations have 
an important influence. The easier way of estimating it 
is via extrapolation (in the rank-plot) of a linear fit to 
the logarithm of the lowest populations. A more sophis- 
ticated method employs a non-linear fitting of the raw 
data to the MaxEnt rank-distribution via this single pa- 
rameter, using Eq. ( 10 1 together with the definition of 
A Eq. (j8J). Such simple procedure in the log-scale works 
nicely enough. Note that log(xo) is just a shift of the dis- 
tribution as a whole, playing no role in its actual shape. 
We have obtained log(xo) = 4.826 (xq = 125), with an 
standard error of 0.036 and a correlation coefficient of 
R = 0.99968. The comparison of the raw data with the 
analytical rank-distribution is displayed in the top panel 
of Fig. [I] The shadowed area represents the 90% confi- 
dence interval due to finite size effects for n c = 141, nu- 
merically estimated, as described in the Appendix. Since 
the empirical rank-distribution falls inside the confidence 
interval, and the observed dynamics are compatible with 
the dynamical assumption of Eq. Q, we consider that 
this case constitutes strong evidence for the applicability 
of the MaxEnt principle. 



2. More detailed analysis 

We have applied the same (zeroth order) approximation 
of proportional growth to the rest of Spain's provinces, 
also including the values of the total populations N in 
the fitting procedure. Using the corresponding values 
of n c together with the estimated values of xq and N 
we get the value for A, and apply then the Gamma- 
scaling law described in the previous Section. The scaled 
rank-distributions are displayed in Fig. 2 (raw data in 
the inset), and the numerical values of the fitted re- 
sults can be found in the Additional Material. A 
quite nice adjustment ensues in general (correlation co- 
efficient of R = 0.99984 for Tenerife) but not in all in- 
stances (R — 0.99635 for Guadalajara). This "failure" 
is linked to the strong correlations between ii and u for 
some provinces, that reaches significant values (of up to 
R = 0.57 for the province of Lugo). Although such cor- 
relations may compromise the validity of the zeroth or- 
der approximation of proportional growth, the quite nice 
scaling-features exhibited by this plot are remarkable in- 
deed. 

It is worth mentioning that the fitted value for the total 
population is found to be systematically lower than the 
actual value for it. Such scenario usually changes when 
the capital city of the province is not considered in the 
fitting process, and thus in the estimation of N . This 
indicates that the population of the capital city is sys- 
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FIG. 1: Top panel: Rank-distribution of the province of 
Alicante (2010), compared with the MaxEnt prediction for 
proportional growth using the actual total population N = 
1926285 and the estimation xo = 125. The shadowed area 
represents the 90% confidence interval for n c = 141 (2010). 
Inset: relative change vs. natural logarithm of the popula- 
tion for this province for the last 15 years, where each gray 
dot represents one municipality for each annual period. Black 
dots denote the mean value in each interval Am = 1 and the 
error bars the standard deviation at that interval. Bottom 
panel: Gamma-scaled rank-distribution of all 50 provinces. 
Inset: raw data (2010). 



parameter log(x ) = 10.80 ± 0.05 (R = 0.99989, higher 
than that of Alicante) and encounter that all the 50 cities 
get located inside the confidence level. When studying 
the relative increment u vs. the log-size u, we find a very 
low size-dependencies (R = 0.15, lower than the one pre- 
vailing for Alicante). This observation confirms the ap- 
propriateness of a proportional growth dynamics. Again, 
the capital city of Madrid exhibits a larger-than-expected 
population (just in the limit of the 90% confidence level), 
indicating the possible existence of a higher, international 
system of correlations. 



C. Higher order approximation 

1. Preliminaries 

Let us proceed now to tackle a power-law growth by as- 
suming some kind of dependence of (u(u)) on u. Starting 
from the simplest, linear relation, the idea is to consider 
an exponential function and its expansion 



bu ~ k\ 



k q e 



(q-l)l. 



(26) 



which is equivalent to Eq. (22 1. We thus face a "q- 
dynamics" with proportional noise, whose distribution 



derives from MaxEnt for q ^ 1 Eq. ( 23 ) 



We confirm, a posteriori, that this procedure is adequate 
up to q = 2, evidencing more information than what 
one would get from just a first order relationship. The 
actual q— value can be determined directly either from 
the dynamical data or from the rank-distributions. We 
can demonstrate via MaxEnt that a cause-effect relation- 
ship exists between the dynamics and the city-population 
distribution, if the two kinds of estimations do match 
(within standard errors). 

We have estimated the parameters xq, A, q, and a of 
Eq. ( |23[ ) for each province from the rank-distribution 
using the methods described in the Appendix, and also, 
now from the dynamics, the q- value for the last 15 years. 
We find (mean value ± standard deviation): log(xo) ~ 

5.5 ± 1.3; log(A) 4.3 ± 1.6; q ~ 1.20 ± 0.45; and 

a ~ 0.43 ±0.24. 



tematically larger than what one would expect from the 
MaxEnt prediction (off the 90% confidence level). This 
is not surprising since, as mentioned above, provinces 
are not ideal, isolated systems: the actual administrative 
municipality for capital cities are usually the sum of the 
historical ones plus some near neighbors, and their econ- 
omy is expected to be highly correlated with that of other 
capital cities. Although only 50 cities is a small sample, 
we have studied the rank-distribution of these capitals 
and their dynamics to shed some light onto this observa- 
tion. We have found the pleasing result that these capi- 
tals form a scale-free system of their own. Comparing the 
rank-distribution with the MaxEnt prediction that uses 
i) the actual value N — 15528025 with ii) the only fitting 



2. Casuistics 

As a general trend, we can fit the rank-distribution with 
no too many technical complications. In many cases the 
capital city has to be excluded, as found in the previous 
subsection. Also, in some other instances, a few of the 
largest cities have to be excluded. We have also found 
'outsiders' for very low-population centers. More details 
are described in the Appendix. 

We encounter a strong correlation between the number of 
municipalities n c and the minimum size xq (R = —0.77 
for log(n c ) vs. log(a;o)). In similar fashion as for USA 
(from East to West), the area of the Spanish' territories 
grows form North to South, due to historical reasons. 
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FIG. 2: Rank-distribution for Valladolid, Castellon, Girona, Burgos, Guipuzcoa, Navarra, Leon, Avila, and Barcelona (year 
2010), compared with the pertinent MaxEnt prediction. The shadowed area represents the 90% confidence interval for each n c 
value. Insets: relative growth vs. logarithm of the population, same as inset of top panel in Fig. 1 but with Au — 0.25. 



Proportional drift noise, scaled by a, has shown to be 
negligible in some cases but quite important in others. 
An low correlation is found with log(a:o) {R = —0.15), 
but it is largely dependent on q and log(A) (R = 0.27 
and R = —0.40). Also, the correlation between q and 
log(A) is rather important (R = —0.37). 
The dynamics' fit presents some preliminary difficulties. 
The main one is high noise for low city-sizes (the cause 
can be found in the last term of Eq. Q, that can not 
be neglected). Since the mean value of this noise is zero, 
we have found it compelling to estimate q via the 'local' 
mean value of u in Eq. (22), expressed as a function of 
u (see Appendix) 



(u)(u) = (k 1 ) + (k q )e^ u . 



(27) 



Remarkably enough, we systematically find here a better 
correlation than for the simple linear fitting. In few cases, 
even if the variance of k q is large, troubles arise when the 
value of (k q ) is small in relation to the drift. In these 
cases the convergence of the fit falls down to the generic 
result q = 1 and an accurate estimation of q cannot be 



achieved from the dynamics. 



3. Results 

We show in Fig. [2] some interesting examples of rank- 
distributions (plots for all provinces are found in the ad- 
ditional material^, and also the table with all the nu- 
merical values). Most of the provinces follow quite well 
the analytical curve, with very few exceptions (the most 
dramatic of those exceptions correspond to Salamanca, 
Orense, and Zamora, where a few cities account for the 
main part of the province's total population and the small 
villages follow a log-normal distribution). In general, bet- 
ter fits are obtained for those provinces for which the 
distributions' evolution during the last 15 years has been 
smooth and slow. High rates of change correlate with 
high a. This rather surprising result tells us that if the 
system is able to reach dynamic equilibrium, it converges 
to the MaxEnt prediction. We plot in Fig[3]the evolution 
of q for some examples, that clearly display the correla- 
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FIG. 3: Evolution of q via MaxEnt distributions for Ali- 
cante (circles), Guipuzcoa (squares), Girona (diamonds), and 
Navarra (triangles). Alicante is the province with larger a 
(0.15) of these four examples, correlated with its larger time- 
variation. 
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FIG. 4: Comparison of measures of q, via the comparison 
MaxEnt distribution vs. dynamics, (gray color represents 
those cases where the dynamical fit fails, see text). 



IV. FURTHER REMARKS 

A. Origin of the g-dynamics 

In the previous section we have successfully parameter- 
ized the dynamics without discussing possible underlying 
mechanisms. We have encountered two general situa- 
tions: 

• i) migration from small villages to big cities and 

• ii) saturation of big cities. 

The cases of Avila, Castellon, Valladolid and Navarra 
constitute nice examples of the first scenario, the most 
frequent for Spain's provinces. We usually find for them 
a value q > 1. Many examples share for few of the 
largest cities the second possibility, with a relative growth 
lower than expected. An explanation for this situation 
can be found in the finiteness of the resources to make 
grow a city, thusly avoiding arbitrary large growth rates. 
Anyhow we have found some examples, as shown for 
Guipuzcoa, where all the cities follow the second trend 
and finding values of q < 1. These cases clearly indicate 
a generalized migration form cities to small towns, due 
to some local socioeconomic paradigm. 

In other remarkable cases both situations arise in the 
same region, as happens for Madrid and Barcelona. This 
is a mark of non-monotonic behavior for u with respect 
u, which in such instance compromises the dynamical fit 
(shown here for Barcelona and in the Additional Material 
for MadriopS). On the other hand, the case of Girona 
deserves some words. Even if it presents a tiny migration- 
tendency, the system is very near total equilibrium, with 
an excellent fit for the rank-distribution. 
The use of power laws to model both migration and satu- 
ration is here of a heuristic nature. We have no evidence 
of any 'microscopic' mechanism generating them, and the 
value of q is merely obtained from empirical observation. 
Identifying the mechanism that generates the dynamics 
constitutes a formidable challenge. It is worth mention- 
ing that growth with saturation has been traditionally 
modelled using the so-called logistic function, [or its gen- 
eralization Y(t)] with parameters a, K, and v that follow 
a differential equation of the type^l 



Y{t) = a{l-(Y{t)/KY}Y{t), 



(28) 



tion we are concerned with here. 

We depict in Fig. [4] the comparison between both in- 
dependent sources for g, MaxEnt vs. dynamics, which 
represent the main result of this work (we display in grey 
the cases where the dynamical fit fails, near the dynam- 
ical value q = 1). We find a mean proportionality of 
0.85 with a correlation of R — 0.97. In such instances 
the dynamical process determines the rank- distribution of 
the city-population. The equilibrium distribution is that 
found via maximization of the Shannon entropy in the 
terms discussed above. 



which formally is identical to Eq. (22). Accordingly, 



migration and growth with saturation generates hyper- 
cxponcntial growth, but, as far as we know, a pure theo- 
retical determination of the q-value from the underlying 
mechanism is still unknown. 



B. The equation of state 

We have shown that the quantities q, A, and a exhibit 
strong correlations amongst themselves, but not with xq. 
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FIG. 5: Same as Fig. 2 for the state of Ohio (USA) in 2010 
(top), and the general electoral results in United Kingdom in 
2005 (bottom). 



MaxEnt gives a relationship between these parameter in 
the form of Eq. ( 24 ) , which might be read as an equation 



of state that relates all of them. The lack of an important 
correlation with xq is explained because it is a propor- 
tionally constant, being just a shift in the log-scale. Since 
we are dealing with scale invariant systems, the physics 
does not depend on scale, and the equilibrium distribu- 
tions respect this fact. The above correlations may indi- 
cate the appropriateness of using such equation of state 
in future city-population research and thermodynamic- 
oriented analysis. 



terns, as electoral results, in the present and past years. 
As an example we plot in the top panel of Fig. [5] the 
rank-distribution with theoretical fit for the state of Ohio 
(log(A) = -3.7, q = 1.62, log(x ) = 10.1 and a = 0.36). 
We also show in the bottom panel the electoral results for 
the 2005 UK general elections (log(A) = -13.5, q = 1.30, 
log(x ) = 4.15 and a = 0.05). Deeper understanding of 
each territory /election, correlated with the kind of study 
presented here may help to better appreciate the effects 
of regional policies. Work along such lines is in progress. 



V. SUMMARY 

We have shown that the observed city-population dis- 
tributions of the Spanish provinces follow in general the 
predicted MaxEnt equilibrium distributions, according to 
their intrinsic growth dynamics. 

• We first have considered a zero-th order approach 
to the problem assuming proportional growth, thus 
finding that the empirical distributions can be 
nicely scaled following the Gamma Scaling Law, de- 
rived from the equilibrium distribution of an scale- 
free system. 

• Secondly, we have considered hyper-exponential 
growth with some proportional drift. We can ac- 
count in this way for both migration and 'natural' 
growth of the population, obtaining better fits than 
in the first case. 

We have checked that the value of the exponent for 
the hyper-exponential growth, estimated from the rank- 
distributions, is equivalent to the one estimated directly 
from the dynamics, which we read as a confirmation of 
the validity of our approach. 



C. Beyond Spanish cities' population 

We have also found evidences of MaxEnt equilibrium 
both in other countries and in other kind of social sys- 
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Appendix A: Estimation of the MaxEnt distribution parameters from empirical rank-distributions 
1. Estimation of the distribution's 90% confidence interval 

The confidence levels of the rank-distributions for a given number n c are estimated as follows: 

i) A list of n c random numbers is generated, following the desired distribution, by inverse transform sampling. 

ii) The list is sorted from largest to lowest n c in order to obtain the rank-distribution. The list is saved for further 

use. 

iii) A large number of lists is generated following i) and ii), obtaining a distribution of numbers for each element of 

the list. 

iv) The 0.95-th and 0.05-th quartiles are obtained from the distribution for each element, determining the lower and 

upper limits of the 90% confidence interval. 

2. Fit for q-exponential growth 

A first method to estimate the MaxEn t q- distributions' parameters A, xq, and q for a given empirical data is a 



direct fit of the rank-distribution to Eq. (20 1. We have used the Mathematica software®! to this end, by means of 
the NonlinearModelFit function using different guesses for the initial values. We also compare the pertinent results 
with a solution obtained from the reproduction of the first three moments of the distribution, finding, statistically and 
numerically, better stability using the logarithmic moments ((\og(x/xo)) n ) instead of ((x/xo) n ). We have tabulated 
them as functions of q and A in the form M n (q, A) = ([log(x/xo)]™), which are defined via 



(log(x/x o) ) = w — - Gl o 
([log(x/x )] 2 > = — 2 - Gil 



1 


r(i 


-q, A) 




2 


r(i 


A) 




n! 


r(i 


-g,A) 



1 1 

(1 - q) 
1 1 1 

(1 -?) / (Al) 



1 ... 1 

... (l- 9 ), 



where G stands for the so-called Meijer G-functions. The associated system of equations is 

M n (q,A) = £[log"(x/x )] - ]T (-l) m ( I, ) E[\og n - m (x)} log m (x ), (A2) 

m=0 ^ ' 

that we deterministically solve using n = 1,2,3 via A, Xo, and q, with the empirical expected values £'[log n (x)] = 

Er=i io s n (^)/« c - 

We have also introduced into our equations the proportional drift parameterized with a. Its inclusion in the above 
equations can be easily materialized using the additive property of cumulants for convolutions. Remember that the 
cumulants of a probability distribution (PD) are a set of quantities that provide an alternative to the PD-moments. 
These, in turn, determine the cumulants in the sense that any two probability distributions whose moments are 
identical will have identical cumulants as well, and similarly the cumulants determine the moments. In some cases 
theoretical treatments of problems in terms of cumulants are simpler than those using moments.^ Taking into account 
that for the centered normal distribution only the second cumulant does not vanish, we face the following system of 
equations 

^7> 2 *M„_ 2l (,z,A) = ]T(-l) m n m £[log"- m (x)]log m (x ), (A3) 
i=0 m=0 ^ ' 

where T™ is the i-th element of the n-th row of the triangle of Bessel numbers® T™ = n\/(i&{n — 2i)\). 
The ensuing system of equations is solved using the FindRoot functionality of Mathematica with different starting 
points so as to ensure the best possible result. We accept a set of parameters if the three different results (fit, 
and equations' system with/without drift) are reasonably similar. When that does not happen, we look for outsiders. 
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excluding some points from the tails (largest or smallest) until reaching a satisfactory convergency of the three results. 
The most frequently outsiders-case obtains for the capital city. The second case refers to a few small villages with 
undersized population. All outsiders are indicated in the Additional Material, together with the results of the three 
estimations. On the other hand, as commented in the text, we have found three notorious cases -Salamanca, Orense, 
and Zamora- where no satisfactory outcome has been achieved with this procedure. 



3. Fitting the dynamical data 



As commented in the text, we have fitted the mean value (u){u) to Eq. (27) via (fci), (k q ), and q. We have employed 
again the NonlinearModelFit function, with different guesses for the initial values. The data set for (u)(u) (and 
also the standard deviation) is systematically computed in bins of size Am = 0.25 for all the provinces. Very few 
points can be found in some of the bins, which introduces high numerical error. Accordingly, we include the bin in 
the data set if, and only if, a minimal number of points exists. We have assumed, for all the provinces, that this 
number is the 15% of the bin with the larger number of points. The bins used in each fit are shown in the inset of 
the rank-distribution plots. Anyhow, this filter is not enough in a few cases, for which a satisfactory result can not 
be encountered. A q-value closed to 1, albeit ill-defined, is the result. 
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